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Abstract 

Consider being given a mapping ip : 

S d-i L4 90, with dQ the 

onto 

(d — l)-dimensional smooth boundary surface for a bounded open simply- 
connected region Q, in R d , d > 2. We consider the problem of constructing 
an extension $ : Bd — > fi with Bd the open unit ball in R d . The map- 

onto 

ping is also required to be continuously differentiable with a non-singular 
Jacobian matrix at all points. We discuss ways of obtaining initial guesses 
for such a mapping $ and of then improving it by an iteration method. 



1 Introduction 

Consider the following problem. We are given 

ip : dB d ^4 dn (1) 

onto 

Notationally, d > 2, B d is the open unit ball in R d with boundary = 
dBd, and is an open, bounded, simply-connected region in M d . We want to 
construct a continuously differentiable extension 

$ : B d i4 H (2) 

onto 



such that 



*ls««-i=P (3) 
det(D$(x)) ^ 0, xeB d (4) 



£><!> denotes the dxd Jacobian of <£>, 



As a particular case, let d = 2 and consider extending a smooth mapping 



ip ■. s 1 i-4 on 

onto 

with fl an open, bounded region in R 2 and (p a smooth mappinng. In this case, 
a conformal mapping will give a desirable solution $; but finding the conformal 
mapping is often nontrivial.. In addition, our eventual applications need the 
Jacobian Z?$ (see [T], [3], [5]), and obtaining _D$ is difficult with most methods 
for constructing conformal mappings. As an example, let tp define an ellipse, 

<p(cos6»,sin6») = (a cos 9, b sin 9) , < 9 < 2ir 

with a, b > 0. The conformal mapping has a complicated construction requiring 
elliptic functions (e.g. see [21 §5]), whereas the much simpler mapping 

$ (xi,x 2 ) = (axi,bx 2 ) , x S B 2 

is sufficient for most applications. Also, for d > 2, constructing a conformal 
mapping is no longer an option. 

In §2 we consider various methods that can be used to construct $, with 
much of our work considering regions fl that are 'star-like' with respect to the 
origin: 

ip(x) = p(x)x, xeS" (5) 
p : S*- 1 ±4 M> 

onto 

For convex regions O, an integration based formula is given, analyzed, and 
illustrated in §3. In §4 we present an optimization based iteration method for 
improving 'initial guesses' for $. Most of the presentation will be for the planar 
case (d = 2); the case of d = 3 is presented in 



2 Constructions of $ 

Let O be star-like with respect to the origin. We begin with an illustration of 
an apparently simple construction that does not work in most cases. Consider 
that our initial mapping ip is of the form (|5|) . Define 

®(x,y) =rp(9) (cos 0, sin 0), 0<r<l (6) 
= p{9)x (7) 

with x = (r cos 0, r sin 9) and p(0) — p(x) a periodic nonzero positive function 
over [0,27r]. This mapping $ has differentiability problems at the origin (0,0). 
To see this, we need to find the derivatives of p(6) with respect to x and y. Use 

9 = tan -1 (y/x) , y > 0, x ^ 



2 



and an appropriate modification for points (x, y) in the lower half-plane. We 
find the derivatives of p using 

dp(0) ^ 86 dp{6) ae 



Then 



dx dx ' dy dy 

86 -y 86 x 



Using these, 



dx x 2 + y 2 ' dy x 2 + y 2 



dx \ x 2 + y 2 ' x 2 + y 2 

1 1 ^ { 6),p{6) + -^-,p/{6) 



The functions 



dy \x 2 + y 2 ' x 2 + y 2 



9 9 

x z xy y z 



x 2 _|_ y2 ' x 2 _|_ yi ' x 2 + y 2 

are not continuous at the origin. This concludes our demonstration that the 
extension $ of dU) does not work. 



2.1 Harmonic mappings 

As our first construction method for $, consider the more general problem 
of extending to all of Bd a real or complex valued function / defined on the 
boundary of Bd- Expand / in a Fourier series, 

n=l 

Define F on B using 

F (x) = -a + r n [a„ cos (n0) + 6„ sin (n9)} (9) 

with x = (r cos 0, r sin 0) . Note that this is the solution to the Dirichlet problem 
for Laplace's equation on the unit disk, with the boundary data given by /(0), 
< < 2tt. 

It is straightforward to show that F is infinitely differentiable for |x| < 1, a 
well-known result. In particular, 

8F °° 

— — = a\ + > (m + 1) r m [a m+ i costo0 + b m+ i smmd] (10) 
dx\ ' 

m— 1 

<9F °° 

t; — = &i + ^ (m + 1) r m [— a OT +i sinm0 + b m+1 costo0] (11) 
ox 2 ^ 
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Figure 1: Starlike region with p denned by (1131) with a = 5 



Depending on the speed of convergence of ((U), we have the partial derivatives 
of F(x) are continuous over B2- In particular, if we have 



n=l 



n \a„\ < oo, 



n=l 



n \b n \ < oo 



then dF/dxi and dF/dx2 are continuous over £?2- 
Given a boundary function 



<p(6) = (<Pi (0),<P2 (*)) 



< 9 < 2tt, 



(12) 



we can expand each component to all of B2 using the above construction in (|9]), 
obtaining a function $ defined on Bd into R 2 . A similar construction can be 
used for higher dimensions using an expansion with spherical harmonics. It is 
unknown whether the mapping $ obtained in this way is a one-to-one mapping 
from B2 onto fi, even if fi is convex. 

The method can be implemented as follows. 

• Truncate the Fourier series for each of the functions ipk(9), k — 1,2, say 
to trigonometric polynomials of degree n. 

• Approximate the Fourier coefficients {a,j} and {bj} for the truncated se- 
ries. 

• Define the extensions (x) in analogy with ^ . 
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Figure 2: The Jacobian for JTSJ with 0.905 < |det (D$ (x))\ < 75.314 



Example 1 Choose 

p{9) = a + cos 6 + 2 sin 29 (13) 

with a chosen greater than the maximum of \cos8 + 2sin2f?| for < 8 < 2tt, 
approximately 2.2361. Note that p(9)cos6 and p(9)sm9 are trig polynomials 
of degree 3. Begin by choosing a — 5. Choosing n — 3, we obtain the graphs 
in Figures\]\ and\^ Figure\]\ demonstrates the mapping by showing the images 
in Q of the circles r = j/p, j = 0, , . . . ,p and the azimuthal lines 9 — irj/p, 
j = 1, . . . , 2p, p = 15. For the numerical evaluation of the Fourier coefficients, 
the trapezoidal rule with 10 nodes was used. Figured shows |det (D$ (x))\ The 
figures illustrate that this $ is a satisfactory mapping. However, it is possible 
to improve on this mapping in the sense of reducing the ratio 

max |det (£><!> (.t))| 

A ($) = ^2 (14) 

1 ' ~ mm |det(D$(x))| 1 ' 

xeB 2 

For the present case, A ($) = 100.7. An iteration method for decreasing the size 
of A ($) is discussed in As a side-note, in the planar graphics throughout 
this paper we label the axes over the unit disk as x and y, and over Q, we label 
them as s and t. 

In contrast to this example, when choosing a — 3 in (I13[) the mapping $ 
derived in the same manner is neither one-to-one nor onto. Another method is 
needed to generate a mapping $ which satisfies (fT3|) J2])-(|4]). 
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Figure 3: Star like Cassini region with p denned in <\11\ with a = 1.5 



2.2 Using C°°-modification functions 

Let (x,y) = (r cos 8, r sin 6), < r < 1. As earlier in (O, consider f2 as star-like 
with respect to the origin Introduce the function 

r(r;«) =exp (k (l - ^j) , 0<r<l (15) 

with k > and T(0, k) = 0. Define $ by 

a = $(x;fc,w) = [T(r;«)jo(0) + (l-T(r;re))w]x, x e B 2 (16) 

with a; = r (cos#, sin6*), for < r < 1, < 6* < 2-7T, with some cj > 0. This is an 
attempt to fix the lack of differentiability at (0,0) of the mapping ©-(O. As r 
decreases to 0, we have $ (x) w lux. Thus the Jacobian of $ is nonzero around 
(0,0). The constants k, oj are to be used as additional design parameters. 

The number uj should be chosen so as to also ensure the mapping $ is 1-1 
and onto. Begin by finding a disk centered at (0, 0) that is entirely included in 
the open set f2, and say its radius is cj , or define 

ujn = min p(6) 

O<0<2tt 

Then choose u> € (0, ljq). To see this is satisfactory, write 

$ (x; k, id) = f(r) (cos 8, sin 8) 
f(r) = r [T (r; k) p{8) + (1 - T (r; «)) u] 



G 



fixing 9 <E [0,27r]. Immediately, / (0) = 0, / (1) = p{9). By a straightforward 
computation, 

/'(r) = -{T[(jfr- W )(r + l)]+ W } 
r 

where T = T (r; k) and p = p (0). The assumption < w < ujq then implies 

/' (r) > 0, < r < 1 

Thus the mapping / : [0, 1] — > [0, p(9)] is 1-1 and onto, and from this $ : B2 — >• £1 
is 1-1 and onto for the definition in (fT6)) . 

This dehnition of $ satisfies ([2])-(j4|), but often leads to a large value for the 
ratio A ($) of (p~4]l . It can be used as an initial choice for a $ that can be 
improved by the iteration method defined in §4. 

Example 2 Consider the starlike region with 

p{9) = \jcos{29) + ^Ja-sm 1 (20) (17) 

with a > 1. The region fi is called an 'ovals of Cassini'. We give an example 
with a = 1.5, (k, uj) — (1.0,0.5). Figure [5| is the analog of Figured For the 
Jacobian, 

minD$(x,y) = 0.0625 

r<l 

max D$(x,y) = 4.0766 

r<l 

The ratio A ($) = 65.2 is large and can be made smaller; see Example \14\ 

A variation to (fTo) begins by finding a closed disk about the origin that 
is contained wholly in the interior of Q. Say the closed disk is of radius 5, 
< S < 1. Then define 

$ (x; k, uj) 

x, < r < S , N 

(18) 

where x = r (cos 9, s'm9). Then the Jacobian D$ around the origin is simply 
the identity matrix, and this ensures that detD$(x) ^ for x 6 E>d- Experi- 
mentation is recommended on the use of either (fl6]) or (fT8]). including how to 
choose k, ui, and S. 

The methods of this section generalize easily to the determination of an 
extension $ : B3 - — \ Q for the given boundary mapping 

onto 



onto 



Examples of such are illustrated later in 
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3 An integration-based mapping formula 



Begin by considering a point P = r(cosa,sina) G P2, r G [0,1), a G [0, 2%). 
Given an angle a < 9 < tt + a, draw a line L through P at an angle of 9 with 
respect to the positive xi-axis. Let P+{9) and P-(9) denote the intersection of 
this line with the unit disk. These points will have the form 



P + (9) = P + r + (0) V , 
P- (9) = P- r_ (0) T]. 



(19) 



with 

r] = (cos 9, sin 9) , a < 9 < tt + a. (20) 
We choose r + (9) and r_ (#) to be such that 

\P+(9)\ = \P + r + (9) V \ = 1, \P-(6)\ = \P - r_ (0) r,| = 1 

and 

|r + (<?)| = |P-P + ((9)|, |r_(0)| = |P-P_(0)|. 

Define 

(21) 

using linear interpolation along the line L. Here and in the following we always 
identify the function (p on the boundary of the unit disk with a 2tt periodic 
function on the real number line. Then define 

$ (P) = - f + (9) d9 (22) 

J a 

We study the construction and properties of $ in the following two sections. 



3.1 Constructing $ 

The most important construction is the calculation of P+ (9) and P_ (9). We 
want to find two points 7 that are the intersection of dP>i and the straight line 
L through P in the direction r], \tj\ = 1. Since P £ int (P2), we have \P\ < 1. 
We want to find 

-y = P + rq, [*y] = 1 

with 77 denoting the direction from P as noted earlier. With the assumption 
(|20|) on 77, we have 

0<P-77<|P|, a<6»<a+i?r 
0<-P-77<|P|, a+i7r<6l<a + 7r 

Using 7-7 = 1, 

\P + rrj\ 2 = P-P + 2rP ■ 77 + r 2 = 1 
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r 2 + 2rP T] + P ■ P - 1 = 

<o 

which implies the roots are real and nonzero. Thus the formula 



f) ■ P± \I(P ■ T]) 2 + 1 - P ■ P 



= r cos(0 — a) ± y 1 — r 2 sin 2 ( 

defines two real roots. Here we see that 

\P- V \< \P\ 
(P ■ rif - P ■ P < 

(P ■ T] f + 1 - P ■ P < 1 

So 



= P-r)+\l(P-r 1 f + l-P-P 



i) + \Jl-r 2 sin 2 (6>- a) 



= -P ■ rj + y (P ■ f]f + 1 - P ■ P 

= —r cos(# — a) + \Jl — r 2 sm 2 (9 - a) 



r_ + r + = 2 \l (P ■ r]) 2 + 1 - P ■ P 



2\ 1-r 2 sin 2 ( 



It is immediate that 



a) 

and therefore the denominator in the formula (|21|) for tp* (9) is zero if and only 
if \P\ = 1 and P _L 77, a case not allowed in our construction. 

Using r_ and r + in (fT9|) . we can construct tp* (9) using (|2"T|) . and this is then 
used in obtaining the mapping $ (P) of (|2"2")l . This formula is approximated 
using numerical integration with the trapezoidal rule. We illustrate this later 
in the section. 

To further simplify the analysis of the mapping $ of (|2"21 , we assume for a 
moment that a — 0, so the point P is located on the positive x-axis. Our next 
goal is to determine the respective angles between P+{9) and P-{6) and the 
positive x-axis. We denote these angles by ip + and ip~, respectively. Using the 
law of cosines in the triangle given by the origin, P, and P + we obtain 

(r + ) 2 = r 2 + 1 - 2rcos(> + ) 
2rcos(> + ) = r 2 + 1 - (r+) 2 

= r 2 + 1 - (-rcos(6>) + y/l - r 2 sin 2 



= 2r 2 sin 2 9 + 2r cos(0) VT - r 2 sin 2 6» 
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- 



-0.5 ■ 




-1 -0.5 0.5 1 1.5 2 



Figure 4: The mapping $ for boundary (f29|) with a = 0.9 



which implies 

V-'+ = i>+(r,< 



arccos (r sin 2 9 + cos(#)\/ 1 — r 2 sin 2 9j 
where we use the function arccos : [—1,1] >-> [0,7r]. Similarly we get 



(23) 



V»_ = r/ 7 - ( r 7 #) = arccos f r sin 2 - cos(0) yl - r 2 sin 2 9 J (24) 

where we use the function arccos : [—1,1] i-> [7r, 27r], 

arccos(ir) = 27r — arccos(a;) 

Using the functions i/>_ and -0+ we can rewrite tp*(6), see pip, in the following 
way: 



1 



<P*{r,9) = -\ I 
1 



■ cos(#) 



2 \ yl — r 2 sin 2 1 



1 - 



■cos(0) 



\/l — r 2 sin 2 i 
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Figure 5: The mapping $ for boundary (|1 T[) with a = 1.5 



This allows us to write formula (|22[) more explicitly in the following way: 



27rJ \ \/l — r 2 s 



sin 



2n Jo V Vi - r 2 sin 2 6» / 



r cos(6>) 



2"" Jo V Vl — r- sni 



2tt 



1 f I rcos(9~ tt) 

2rr/ - ' 7 l-r 2 sin 2 f6»-7r) 



V {i) + {r,6))d6 



i rl rcosje) \ { ^ e))d9 

2lr Jo I Vl-r 2 sin 2 0/ 



If 271 / r cos(d) 



27T Jtv \ \Jl-r 2 sin 2 1 



IT (1+ / C ° S(g) . W*M))^ (25) 



2?r "'o \ Vl-^sin 2 ^/ 
Here we used the variable transformation i— > tt + for the second equality and 
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the new definition 



Mr, 8) : 



ip+(r,0), 0<9<tt 
ip-(r,0-n), T7<e<2ir 



(26) 



where the functions ^_ and ip+ are defined in ([2"B"| and (pM)) . We remark, that 
the function -0* : [0, 1) x [0, 2tt] [0, 2ir] is a continuous function which follows 
from its geometric construction. We further define 



a 27T periodic continuous function on [0, 1) x [0, 2tt]. If we now go back to the 
general case P — r(cos(a), sin(a)), a £ [0, 2w), we can rotate the given boundary 
function ip and obtain 



Before we study the properties of the extension operator £, we present two 
numerical examples. 

To obtain $ (P) , we apply the trapezoidal rule to approximate the integral 
in (|28| or (|22|) . The number of integration nodes should be chosen sufficiently 
large, although experimentation is needed to determine an adequate choice. 

Example 3 Consider 

tp {cos 9, sin 9) = (cos 9 -sin 9 + a cos 2 9, cos 9 + sin 9) , < 9 < 2tt (29) 

with < a < 1. We choose a = 0.9 and apply the above with n = 100 subdivi- 
sions for the trapezoidal rule to evaluate Figure [7] shows the mapping <I>, 
done in the same manner as earlier with Figures^ and\3^ 

Example 4 We consider again the ovals of Cassini region with boundary given 
in Jj?| ) with a = 1.5. The mapping A22\) is illustrated in Figured However, 
for a somewhat closer to 1, this integration formula \22\) no longer produces a 
satisfactory $. 

3.2 Properties of Sip 

To study the properties of the extension operator £, see (|28|) . we have to study 
the behavior of the functions and k, see (12^1) and (f2~T)) . at the boundary r = 1. 
We start with the function and define the values of this function for r = 1 
first: 



k{r,9) := 1 + 



r cos(6*) 



(27) 



\/l — r 2 sin 2 9 




(28) 




(30) 
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Because of 



lim 1 — r = 

i — >i 



and the boundedness of the sine function, the limit 

lim 1 - r 2 sin 2 (6>) = 1 - sin 2 (0), 

i — >i 

is uniform for 9 € [0, ^7r] . The uniform continuity of the square root function 
implies that 



sin 2 i 



lim (cos(0)\/ 1 - r 2 sin 2 o) = cos(0) \J\ 

= cos 2 

uniformly for 9 G [0, ^7r] . Together with similar arguments for the function 
we get 



lim (r sin 2 9 + cos ((9) Vl - r 2 sin 2 9) = sin 2 9 + cos 2 9 = 1 



uniformly in 0. Finally we use the uniform continuity of arccos(-) to conclude 
that 



lim ij)*{r, 0) = 1™ arccos (r sin 2 9 + cos(0) \f\ 

r— yl~ r— ^ 

= arccos(l) 

= =V*(M) 

converges uniformly on [0, ^7r] . Because 



r 2 sin 2 i 



'l -sin 2 (0) = -cob(0), 9e[%ir,%], 
we see in a similar way that 

lim ip*(r, 9)= lim arccos (V sin 2 9 + cos(0)\/ 1 — r 2 sin 2 

= arccos(sin 2 9 + cos(0)(— cos(0))) 
= arccos(— cos(20)) 
= arccos(cos(20 — 7r)) 
= 29 -n 

= V».(M) 

uniformly for 9 € [W, 7rl . Similar arguments apply for € [it, 2tt] and we finally 
conclude that ip*(r,9) converges uniformly to ^>*(1,0) as r approaches 1. This 
proves the next lemma. 

Lemma 5 The function ip*, defined by \26\) and \30\) , is continuous on [0, 1] x 
[0,2tt]. 
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We remark that continuity on a closed interval implies uniform continuity. 
Now we turn to the function k defined in ([27]) . Here we define k for the value 
r = 1 in the following way 

f 2, < 9 < ±tt 
fc(l,0) := < 0, ±tt < 6> < |tt (31) 
[2, §7T < 6» < 2tt 

Obviously fc(l,-) cannot be the uniform limit of k(r, •) as r approaches 1, but 
the following lemma holds. 

Lemma 6 The function k : [0, 1] x [0, 2n] H> [0, 2], defined by and (T57]), is 
bounded; and for every S > 0, the function k(r,9) approaches k(l,9) uniformly 
on Is as r approaches 1. -Here 

/a := [0, 2tt] \ { (±tt - 5, ±tt + <S) U (|tt - 6, f tt + 5) } 

Proof. That k is bounded by 2 follows from 



V 7 ! - r 2 sin 2 > Vl - sin 2 = | cos(0)| 

and the fact that r € [0,1]. The function is uniformly continuous on 

e < z < 1 for every e > 0. From the proof of Lemma [3] we know that 

lira 1 — r 2 sin 2 9 = cos 2 9 

uniformly for <E [0,27r]. Together with the uniform continuity of l/y/z on 
[cos 2 (^7r — (5), 1], this shows 



l im 1 + rcos ^ = 1 . cos ^) 



y/l - r 2 sin 2 9 \cos(9)\ 
uniformly on Is- Remembering 



|cos(0)| 



cos(0), 9 e [0, §7r] U [§7r,27r] 
cos(6>), 6> e [|tt, §7t] 



proves the Lemma. 



Motivated by the properties of tp* and k we now prove a more general result 
for integral operators of the form (f28|) . 



Lemma 7 Let fci,fc 2 : [0,1] x [0,27r] 1— > R &e bounded functions which are con- 
tinuous on [0, 1) x [0, 27r]. Assume there is a finite set E = {6\, . . . ,9 n } such 
that 

lim ki(r,9) = ki(l,e), *=1,2, 
1 — >i- 
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uniformly on Ig := {9 6 [0, 2ir] | \6 — 9j\ > 5, j = 1, ...,n} for every S > 0. 
Then for a periodic continuous function ip : [0, 2n] t— > R the function 



$(r,a):= / fc x (r, 6)<p(k 2 {r, 9) + a 
Jo 



)d0 



is continuous on [0,1] x [0, 2-7r] and 2n periodic in a. 



Remark 8 The above lemma will apply to each component of the function Sip 
defined in H28\) with k\ = k and k 2 — ip* and E = {^7r, §7r}. This shows the 
continuity of Sip. 

Proof. The uniform convergence on Ig, S > arbitrary, shows that •), 
i = 1,2, are piecewise continuous and bounded functions on [0, 2n], so all in- 
tegrals exist. The continuity of a) on [0, 1) x [0, 2tt] follows easily from 
the continuity of the functions ki, i — 1,2. The periodicity follows from the 
periodicity of p and the definition of $. So we only need to show the continuity 
of ^(r, a) on {1} x [0, 2tt] for example at (1,57). Because of the periodicity of 
<i>(r, •) and the property (Sip) (a) = (Sip a )(0), where ip a (9) — ip(a + 9) we only 
need to prove the continuity for one value of a, for example a = ir. We estimate 



|*(r,a)-3(l,7r)| 



< 



kx (r, 9)p(k 2 (r, 9) + a) - k x (1, 9)<p(k 2 (l,9) + tt) d9 



kx{r,e){(p{k 2 {r,e) + a) - <p(k 2 (l, 9)+ir)) d9 



271 



(h(r,9)-h(l,9)Mk 2 (l,9) + ir))d9 

Now we know that ki, k 2 , and ip are bounded functions, for example 

|fti(r, |fci(r, < M, M> 

for all (r, 0) S [0, 1] x [0, 2ir\. So we only have to show 



lim 

r->l 



2;r 



and 



|^(fe 2 (r,0) + a) - ip(k 2 (l, 9)+n)\d9 = 



lim / \h(r,9) - fci(l,0)|d0 = 



(32) 



(33) 



We start with the first limit. Given an e > we choose 8 > small enough such 
that 



2M d9= - 
2 



(34) 

l[0,2n]\l s 

Now we observe that tp is uniformly continuous on R because it is continuous 
and periodic. So there is a uj > such that 



\cp(x) - <p(y)\ < — 



(35) 
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if \x — y\ < oj. We also know that k 2 (r, •) converges uniformly on Ig to k 2 (l, •)) 
so there is a ro € (0, 1), such that 

| fe2M )-W)|<f 

for all r > ro and 9 E Is. If furthermore |a — 7r| < u>/2, we conclude that 

|(fca(r,0) +a) - (MM) + tt)| < \k 2 (r,6) - *i(l,0)| + |a-7r| 

< w 

which by ([55)) implies 

|^(fc 2 (r,0)+a)-^(fe 2 (l,0)+7r)| < ^ (36) 

for all (r, a) G [r , 1] x [tt - w/2, ?r + w/2] and 9 e Is- Combining ((MJ) and (|3"6l 
we can estimate 

/■2-7T 

|<p(fe»(r, 6») + a) - (p(k 2 (l,0) +n)\d9 



< / \<p(k 2 (r,O)+a)-ip(k 2 {l,0)+w)\d0 

J[0,2tv]\Is 

+ [ \ip{k 2 {r,6) + a)-ip{k 2 {l,0)+%)\dd 



< 2Md9+ —d9 

J[0,2-k]\Is Jls 4?r 

< - + 2tt • — = £ 

2 47T 

for all (r, a) £ [r , 1] x [-zr — w/2, 7r + cj/2], which proves (|3"2")l . 

To prove (l33l) we again choose an arbitrary e > and choose S > such that 
(|34|) is true. Now the uniform convergence of ki(r, •) to fci(l, •) on Is proves the 
existence of a r x g (0, 1) such that 

|AiM)-*i(M)l<^: (37) 

for all (r, 0) € [n, 1] x /j. Using (JMJ) and (J37]) we estimate 



\k 1 {r,9)-k l {l,9)\d9 = / \k 1 (r,6)-k 1 (l,6)\d6 
+ / |*i(r,0)-*i(l,0)|d0 

Jls 

I 2Md9+ —d9 



[0,27r]\/ { J/j 4?r 

£ „ e 

< - + 2tt = £ 

- 2 4tt 

for all re [n, 1] . This proves (|33|) . 

Now we state the results about the extension operator £. 
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Theorem 9 Let<p : <9i?2 H> K 2 , be a continuous function, then^(P) = (£ip)(r,a), 
P G B2, see A28\) , is continuous function on Bi o,nd 



®\dB 2 = <f 



(38) 



Proof. In Lemma [5] and Lemma [5] we have shown that the functions k 
and i/j* in (j2"5)l satisfy the assumptions of Lemma [7J So the continuity of <&(P) 
follows from Lemma [7] For P € 9i?2 the polar coordinates of P are given by 
(r, a) = (1, a), a € [0, 2tt], so we get with dHSI) and ((31) 



Corollary 10 Le£ il C K 2 &e a domain with boundary dfl and tp : dB>2 i-> <9J7 

6e a continuous parametrization of the boundary. Then the function £(p, defined 
in (EE), maps B>2 onto fl. 

Proof. Theorem [9] implies that £ ip : B2 H- R 2 is continuous and that 
(£(p)(dB2) = dQ. We assume that the parametrization (p moves along the 
boundary of Q in the positive direction. For y £ we then have 



where deg is the mapping degree; see [TU1 Chapter 12]. But this implies that 
there is at least one x € B2 such that (£ip)(x) — y. ■ 

Theorem 11 Let fl C M 2 be a convex domain with boundary dfl and tp : 3_B 2 > 
<9f2 be a continuous parametrization of the boundary. Then (£tp){B2) C f2. 




— (TT(p(a) + inp(2ir + a)) 

2lT 

<p(a) = p{P) 



because of the 2n periodicity of tp. 



deg(£V,y) = 1 
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Proof. We have to show (£ip)(P) £ for every P £ B^. We use the first 
equation in formula (|25l) 



» ( p) * r i + ^ i r -,, + „,„, 

" \ W 1 _ t2 cin^ 




j=o \ 2^/1 -r 2 sin 2 (6» J ) / 

/ 2^/l-rW(0 i ) 

where (9j = nj/N and we further assumed again that P is on the positive real 
axis to simplify the notation. Here we have used the fact that the integral is the 
limit of Riemann sums. Each term of the sum is a convex combination of two 
elements of f2 and therefore in f2. But the sum itself is a convex combination, 
so the sum is an element of f2. Finally f2 is closed, so $(P) £ fi. ■ 

The two last results imply that for a convex domain f2 we get (£</?) (-B2) = 
Q, but there is still the possibility that £ (<p) is not injective. Our numerical 
examples seem to indicate that the function is injective for convex fi, but we 
have no proof. For non-convex regions, it works for some but not others. It is 
another option in a toolkit of methods for producing the mapping $. 

The integration-based formula (|2"2"j) can be extended to three dimensions. 
Given 

ip : dB 3 ±4 <9f2, 

onto 

define the interpolation formula (p* (9, uj) as before in (|2ip . with (9, uj) the spher- 
ical coordinates of a direction vector r] through a given point P £ P3. Then 
define 

$(P) = — / / (0,w)sincjdwd0 (39) 

2tt Jo Jo 

A proof of the generalization of Corollary [10] can be given along the same line 
as given above for the disk B^. 



4 Iteration methods 

Some of the methods discussed in §2 lead to a mapping $ in which det (_D$ (x)) 
has a large variation as x ranges over the unit ball Bd, especially those methods 
based on using using the C°°-function T (r, k) of (IT51) . We seek a way to 
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Figure 6: The initial mapping $ for Example [T2l with a = 5, based on (fT6|) 

improve on such a mapping, to obtain a mapping in which det (D<& (x)) has a 
smaller variation over Bd- We continue to look at only the planar problem, while 
keeping in mind the need for a method that generalizes to higher dimensions. 
In this section we introduce an iteration method to produce a mapping <I> with 
each component a multivariate polynomial over B 2 . 

Assume we have an initial guess for our mapping, in the form of a polynomial 
of degree n, 

N„ 

= £a<%(aO, xe B 2 (40) 
j'=i 

We want to calculate an 'improved' value for <J>n°\ call it <£>„. 

The coefficients € M 2 . The functions {tpi, •••,ipN n } are chosen to be 
a basis for II„, the polynomials of degree < n. and we require them to be 
orthonormal with respect to the inner product (•, •) associated with L 2 (B2). 
Note that N n = dim (II n ) = \ (n + 1) (n + 2). As basis functions {ipj} in our 
numerical examples, we use the 'ridge polynomials' of Logan and Shepp [5], an 
easy basis to define and calculate; also see [31 §4.3.1]. 

We use an iterative procedure to seek an approximation 

N„ 

®n{x) = y ^Oi n jll) j {x) (41) 
3=1 

of degree n that is an improvement in some sense on $^ . The degree n used 
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Figure 7: The mapping $ for Example 1121 with a — 5, obtained using iteration 

in defining , and also in defining our improved value , will need to be 
sufficiently large; and usually, n must be determined experimentally. 

The coefficients are normally generated by numerical integration of 

the Fourier coefficients |a^°H, 



„(°) 



where <£> is generated by one of the methods discussed in 
used is 

Zir \ v ^ I 2irm 

1=0 m=0 



g{x,y)dxdy. 
b 2 2P 



(42) 

The quadrature 
(43) 



where g (r, 0) = g (r cos 6, r sin 9). Here the numbers {coi} are the weights of the 
(p + l)-point Gauss-Legendre quadrature formula on [0, 1], and the nodes {n} 
are the corresponding zeros the the degree p + 1 Legendre polynomial on [0, 1]. 
This formula is exact if g is a polynomial of degree < 2p + 1; see [H §2.6]. 

We need to require that our mapping will agree with ip on S 1 , at least 
approximately. To this end, choose a formula q n for the number of points on 
S 1 at which to match $„ with the function cp and then choose {zi, . . . , z qn } on 
S 1 . Require <£>„ to satisfy 



$ n (Zj) = (p (Zj) , j 



(44) 



20 




Figure 8: The boundary for the starlike region with p = 3 + cos 6 + 2 sin 26 



which imposes implicitly q n conditions on the coefficients {o: n ,j}- If y is a 
trigonometric polynomial of degree m, and if n > m with q n = 2n + l, then (|44[) 
will imply that $n|gi — V over dfl. Our numerical examples all use this latter 
choice of q n . 

Next, choose a function J 7 (a), a = [ai, . . . , a/v] , and seek to calculate 
a so as to minimize J- (a) subject to the above constraints (]44[). How should 
J 7 be chosen? To date, the most successful choice experimentally has been 
T (a) — A ($„) , defined earlier in (|T3|) . 

4.1 The iteration algorithm 

Using the constraints (|44|) leads to the system 

= <p, (45) 





V'i(^i) • 








A = 






V = 








■ i>N(z q ) _ 







Because is a trigonometric polynomial of degree n, it is a bad idea to 

have q n > 2n + 1. The maximum row rank of A can be at most 2n + 1. Let 
{zi, . . . , z 9 } denote qvi evenly spaced points on S 1 . We want to minimize T (a) 
subject to the constraints (j45|) . 

We turn our constrained minimization problem into an unconstrained prob- 
lem. Let A = U CV be the singular value decomposition of A; U is an orthogonal 
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Figure 9: The boundary mapping $ for the starlike region 
with p = 3 + cos 8 + 2 sin 28 



matrix of order q, V is an orthogonal matrix of order N, and C is a 'diagonal 
matrix' of order q x N. The constraints (|45|) can be written as 

CVa = U T ip (46) 

Introduce a new variable f3 = Va, or a = V T f3. Then Cf3 = U T ip and we can 
solve explicitly for 7 = . . . , /3 9 ] T . Implicitly this assumes that A has full 

rank. Let S — . . . , /3jv] T , /3 = [7 T , 5 T ] T . Then introduce 

G (S) = T (a) (47) 

using a = V T j3 and the known values of 7. We use our initial |a^°' ) j in (|4T)]) to 
generate the initial value for /3 and thus for S. 

The drawback to this iteration method is the needed storage for the q x N 
matrix A and the matrices produced in its singular value decomposition. In 
the following numerical examples, we minimize G using the Matlab program 
fminunc for unconstrained minimization problems. 

Example 12 Recall Example\]\ with a — 5. Generate an initial mapping $ 
using U6}) with k = .5, uj — 1.0. Next, generate an initial polynomial of 

degree n — 3, using numerical integration of the Fourier coefficients lotj I of 
l{4 2fy . We then use the above iteration method to obtain an improved mapping. 
Figured shows the initial mapping $ 7 and Figure^ shows the final mapping $„ 
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Figure 10: The boundary mapping $ for the starlike region with p from (|17|) 
with a = 1.5 

obtained by the iteration method. With the final mapping, we have 3> n |si = ¥> 
to double precision rounding accuracy, and 

A($) = 6.21 

Compare the latter to A ($) = 100.7 /or the mapping in Example\]\ 

Example 13 Consider again the starlike region using A13\) of Example [IJ but 
now with a = 3. The harmonic mapping of Q2.1\ failed in this case to produce 
a 1-1 mapping. In fact, the boundary is quite ill-behaved in the neighborhood of 
(—0.2,0.2), being almost a corner; see Figure [21 In this case we needed n = 7, 
with this smallest sufficient degree determined experimentally. To generate the 
initial guess $, we used \16\) with (k,lu) — (0.5,0.1). For the initial guess, 

A ($7 0) ) = 840. We iterated first with the Matlab program fminunc. When it 
appeared to converge, we used the resulting minimizer as an initial guess with a 
call to the Matlab routine fminsearch, which is a Nelder-Mead search method. 
When it converged, its minimizer was used again as an initial guess, returning 
to a call on fminunc. Figured shows the final mapping $7 obtained with this 
repeated iteration. For the Jacobian matrix, A ($7) = 177.9, further illustrating 
the ill-behaviour associated with this boundary. As before, = tp to double 

precision rounding accuracy. 

Example 14 Consider again the ovals of Cassini region with boundary given 
in Ji7| ) with a = 1.5. As our initial mapping $, we use the interpolating 
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Figure 11: The optimal mapping $7 for the starlike region with p(0) = 5 
sin 9 + sin 39 — cos 59 



integration-based mapping of 122^ . illustrated in Figured We produce the initial 
guess for the coefficients | a j°''} °f &4*fy by using numerical integration. Unlike 
the preceding three examples, the boundary mapping Lp is not a trigonometric 
polynomial, and thus the interpolating conditions of \4-J$ will not force $ 



""IS 1 



to 



equal if over dfl. For that reason, we use a higher degree than with the preceding 
examples, choosing n = 16. Fiaure \7U shows the resulting mapping $. With this 
mapping, A ($) = 26.11. On the boundary, 

max 1$ (x) - ip (x)\ = 2ME - 4 

xgS 1 

showing the mapping does not depart far from the region f2. 

Example 15 Consider the starlike domain with 

p(9) = 5 + sin0 + sin30 - cos50, < 9 < 2ir. 

in (0|)- ([?}). Using the degree n = 7 and the inital mapping $ based on A16\) 
with (k,lu) = (0.2, 1.4), we obtained the mapping illustrated in Figure [771 The 
minimum value obtained was A ($7) = 6.63. As a side-note of interest, the 
iteration converged to a value of A ($) that varied with the initial choice of (k, uj). 
We have no explanation for this, other than to say that the objective function 
A (<£>) appears to be ill-behaved in some sense that we do not yet understand. 
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4.2 An energy method 

In this section we present a second iteration method, one based on a different 
objective function. Instead of A, see (IT4")) . we use 

x "-' B §§ ii*»(6)-«.a)iis 

Ki Lx 1 

(48) 

We again impose the interpolation conditions given in (|44l) ; and the free param- 
eters are given by 5, see (|47|) . First we explain the definition of the points £j 
and Q appearing in formula (|48|) . The points £i are located inside the unit disk 
and are elements of a rectangular grid 

te|i=l,-A}= (^jnfl 2 ; 

the density of the grid is determined by k\ > 0. The points Q are located on 
the unit circle and distributed uniformly 

{0|i = l,...,ii} = {(c OS (^) ) sin^)) \j = l,..-M 

L\ G N. Furthermore the function A contains the parameter a > 0. So in 
addition to the dimension n of the trial space for $ n , this method uses four 
parameters: q n , the number of interpolation points along the boundary; hi, 
which determines the grid density inside the unit disk; L%, the number of points 
along the boundary; and a, the exponent in formula (1481) . 

The motivation for the function A is the following. We start with an equally 
distributed set of points in the unit disk, {£j | i = 1, . . . ,K±} and we try to 
force the mapping $„ to distribute these points as uniformly as possible in the 
new domain fi. One can think of charged particles which repel each other with 
a certain force. If this force is generated by the potential r~ a then the first 
term in formula (j48[) is proportional to the energy of the charge distribution 
{$n(£t) I i — 1, • • • ,Ki}. When we go back to our original goal of creating 
a mapping $ which is injective, we see that this is included in this fuctional 
because the energy becomes infinite if two particles are moved closer. 

The second goal for our mapping is that $„(i?2) C fi, to enforce this condi- 
tion we use a particle distribution along the boundary of Q given by {& n {Cj) I 
j = 1, . . . , L{\. These charges will repel the charges {$n(C0 | « = 1, ■ ■ ■ , K{\ 
away from the boundary. The energy associated with the interaction between 
the interior points and the boundary points gives us the second term in formula 



So we can consider the algorithm to minimize the function A as an attempt to 
minimize the energy of a particle distribution in Q. This should also guarantee 
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Figure 12: A grid on the unit disk 

that the mapping $„ has a small value for the function A, because the original 
points {£; | i = 1, . . . , K\} are equally distributed. 

In our numerical experiments we used a — 2, so the function A(<£>„) is 
differentiable as a function of the parameter S. Furthermore we adjust k± G N 
in such a way that K\ sa N n and we choose L\ ~ h\. For the parameter q n we 
use the same value as in §4.11 

Example 16 Consider the starlike domain defined in U3\) with a = 5 again. 
We use n = 3, a = 2, K\ — 177 ', L\ — 160. To minimize the junction A we 
use the BFGS method, see jTjj. Figures and \13\ show a rectangular grid in 
the unit disc and its image under the mapping . For the initial guess we 
have A($3°' ) ) rj 11500 and A($3°') sa 29. For the final mapping $3 we obtain 
A($3) 7930 and A($3) w 10. This shows that the function A implicitly also 
minimizes the function A. Figure \T4\ shows the image of the final mapping $. 



5 Mapping in three dimensions 

In this section we describe an algorithm to construct an extension $„ : S3 t— > 17 
for a given function ip : S 2 H> dil. We assume that is starlike with respect 
to the origin. The three dimensional case differs from the algorithm described 
in 2] in several ways. The dimension of n„ of the polynomials of maximal 
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Figure 13: The image of the grid in Figure [T2l under the mapping $ 3 for the 
domain given in (fT3|) . 

degree n is given by N n = (™) , so any optimization algorithm has to deal with 
a larger number of degrees of freedom for a given n when compared to the two 
dimensional case. Whereas in the two dimensional case a plot of «& n (i?2) reveals 
any problems of the constructed $ n with respect to injectivity or $„(i?2) C 
fl a similar plot of <& n (i?3) is not possible. For this reason, at the end of 
each optimization we calculate two measures which help us to decide if the 
constructed is injective and into. 

On the other hand the principal approach to constructing is very similar 
to the algorithm described in JU Again we are looking for a function <i?„ given 
in the following form 

*n0*0 = ^a nJ i>j(x), x £ B 3 , 

3=1 

where {ipi, ■ . ■ , ipN n } is an orthonormal basis of II„ and the vectors a n .j G K 3 , 
j = 1, . . . , N n are determined by an optimization algorithm. 

For a given n £ N we use the extremal points of Womersley, see [9] , on the 
sphere S 2 . We will denote these points by W n = {z[ n \ . . . , z^ +i y}. These 
points guarantee that the smallest singular value of the interpolation matrix 

A n :=: ■ ■ 



27 



Figure 14: The image of the grid in Figure [T2l under the final mapping $3. 



stays above 1 for all n which we have used for our numerical examples. The 
number (n + l) 2 is also the largest possible number of interpolation points on 
the sphere which we can use, because dim(II n |s2) = (n + l) 2 , see [3[ Corollary 
2.20 and formula (2.9)]. Again we enforce 

*„(*<">) = v (z<">), j = l,...,(n+l) a , 
for the mapping function <fr„; see also (j45|) . To define the initial function 

N„ 

*i°)(x) = ^ai°^(x), x & B 3 , (49) 
i=i 

we choose 

«S = j = l...,N n . (50) 

(•, -)b 3 is the usual Li inner product on B 3 . The polynomial 3>^ ' is the orthog- 
onal projection of $ into n„. The function $ is some continuous extension of tp 
to S3 , obtained by the generalization to three dimensions of one of the methods 
discussed in ^2|3I Having determined \ we convert the constrained opti- 
mization of the objective function A(-) into an unconstrained minimization, as 
discussed earlier in (|45"|) -([47 |) . 

Once the Matlab program fminunc returns a local minimum for A («&„) 
and an associated minimizer we need to determine if 4>„ satisfies 

*„(x) 7^ * n (y), x,yeB 3 , x^y, (injective) (51) 
*„(fl 3 ) C Q, (into) (52) 
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For this reason we calculate two measures of our mapping & n . 
Given K £ N we define a grid on the unit sphere, 

S k '■= {( sin (l? 7 ^) cos {-R ilT ) > sin {-R^j) sin (-iwr) ,cos (^ttj)) | 

j = 0,...,K, i = 0,...,2K-l}. 

For L £ N, we define a cubic grid in B3, 

B 3 , L := {j^j n B 3 
so every element in B^ l is given by 

i 2 +j 2 + k 2 < L 2 . 

To measure the minimum of the magnitude of the gradient of ip over S 2 , we 
define an approximation by 

m K (<p) := mm - 

x,yeS 2 K IF -2/11 

This number is used to calculate 

E hK ($ n ) := mm ^ y- / m K (ip) 

x,y£B 3 , L , \\x — y\\ I 

x^y 

Because of * n |s 2 ~ V> we expect Bi,k < 1- We use the occurrence of a very 
small value for £1, to indicate that (|5"Tj) may be violated. The result 
Ei > k(&7i) ~ 1 is the best we can achieve, for example, with ip and 3> n the 
identity mapping. 

If (|52|) is violated there is a point x £ B3 and a point y £ S 2 with 3>„(.t) = 
This shows that the following measure would be close to zero 

2W*n):= -in H $ f-fll U) 

x£B 3}L ,y£S 2 K \\x~y\\ I 

Again we expect E2.k,l{^u) < 1, and a very small value of £?2,K,i(*n) indicates 
that (1521 might be violated. For each <fr„ which we calculate we will always 
report Ei t ic(& n ) and E2,k,l(&ti)- For larger K and L we will get a more 
accurate test of the conditions (|51l) and (|52|) . but the cost of calculation is 
rising, the complexity to calculate £?2,ii",z (*n) for example is 0(n 3 K 2 L 3 ). For 
our numerical results we will use K = 40 and L = 10. 

We consider only starlike examples for f2, with <9f2 given as 

if (x) — p (x) x 1 x £ S 2 

— p{6,4>) (sin 9 cos 4>, sin sin <fr, cos 0) (53) 
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Table 1: Measures of approximation stability for (j5"4")) 



A($ a ) 


£a,4o($e) 


■^2,40,io(^e) 


3.0575574308 


0.7485506872 


0.6626332145 



Table 2: Measures of approximation stability for ([55]) 



Function 


A(-) 




£■2,40,10(0 




394.3717406299 


0.2088413520 


0.5926402745 


$ 6 


43.8782117161 


0.2018029407 


0.5175844592 



p(9,(f>) = p (sin 9 cos <f>, sin sin (/>, cos 9) 

To create an initial guess, we begin with the generalization of ©-([7]) to three 
dimensions, defined in the following way: 

$ (x) = rp(9, <j>) (sin 6 cos </>, sin sin </>, cos 6) 
= p{6,4>)x 

for x = r (sin 6* cos 4>, sin sin <f), cos 0), 0<#<7r, < </> < 27r, < r < 1. 
We assume p : S 2 — > 90 is a given smooth positive function. The initial guess 
is obtained using (|35 ]l -([5D |) . the orthogonal projection of $ into H n . Even 
though $ is not continuously differentiable over B3, its orthogonal projection 
$^°- ) is continuously differentiable, and it turns out to be a suitable initial guess 

with $1 0) rj (p. 

s 2 

Example 17 In our first example we choose 

p(6,4>) :=2 + (cos6>) 2 (54) 

Using n = 6 yields the results given in Table [JJ /or i/ie mapping $g obtained 
using the optimization procedure described above. 

See Figure [731 for an illustration of the images of the various spheres jS 2 . 
In this example the initial mapping <f>^ turned out to be a local optimum, so 
after the first iteration the optimization stopped. The measures E\ and E2 seem 
to indicate that the function is into ft and infective. The error of $g on 
the boundary is zero. 

Example 18 Again the boundary dQ is given by \5S\) . but this time we choose 
p(9, 4>) := 2 + 008 6*+ - sin 9 sin (55) 

Using n = 6 gives us the results shown in Tabled We let $g ^ denote our initial 
guess for the iteration, derived as discussed earlier. 

See Fiaure [Tb] for an illustration of the images of the various spheres jS 2 . In 
this example the A(-) value of the initial mapping is significantly improved 
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by the optimization. During the optimization the measures E\ and Ei do not 
approach zero, which indicates that <E>6 *s o, mapping from B3 into f2 and is 
injective. The error of ^g ' and $6 on the boundary is zero. 
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Figure 15: Images of \S 2 , i = 1, 2, 3,4 
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Figure 16: The images $( 34 ) (\S 2 ) , i= 1,2,3,4 
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